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A  Three-Dimensional  FDTD-PML  Algorithm 
for  Nonlinear  Dispersive  Media 

S.  J.  Yakura 

Air  Force  Research  Laboratory,  Directed  Energy  Directorate 
Kirtland  AFB,  New  Mexico  87117 


Abstract 

Starting  with  the  unsplit-field  uniaxial  PML  formulation,  a  second-order  accurate  FDTD-PML  algorithm  is 
obtained  for  the  first  time  using  the  piecewise-linear  approximation  for  nonlinear  dispersive  PML  media. 

In  the  absence  of  the  PML  interface,  the  nonlinear  dispersive  FDTD-PML  algorithm  reduces  to  the  usual 
nonlinear  dispersive  FDTD  algorithm. 

I.  INTRODUCTION 

With  the  advent  of  high  power  computers  that  provide  fast  execution  times  and  great  quantities  of 
computer  memory,  we  are  at  the  stage  where  we  can  perform  direct  numerical  calculations  of  Maxwell’s 
equations  in  nonlinear  dispersive  materials.  Out  of  many  numerical  techniques  available  in  the 
computational  electromagnetic  community,  one  that  has  shown  a  great  promise  in  the  time  domain  is  the 
well-known  finite-difference  time-domain  (FDTD)  method  [1].  It  is  based  on  using  a  simple  staggered 
differencing  scheme  in  both  time  and  space  to  calculate  the  transient  behavior  of  electromagnetic  field 
quantities.  One  of  the  greatest  challenges  of  the  FDTD  methods  has  been  the  efficient  and  accurate 
formulation  of  electromagnetic  wave  interactions  in  unbounded  regions.  For  such  problems,  an  absorbing 
boundary  condition  must  be  introduced  at  the  outer  layer  boundary  to  simulate  the  extension  of  the  lattice 
to  infinity.  One  approach  that  has  given  a  great  promise  in  realizing  such  an  absorbing  outer  boundary 
inside  the  finite  volume  computational  domain  is  the  well-known  perfectly-matched-layer  (PML)  algorithm 
that  was  first  introduced  by  J.  P.  Berenger  [2]  in  1994  for  the  free  space  boundary.  Since  that  time  Chew 
and  Weedon  [3]  came  up  with  the  modified  PML  algorithm  that  is  based  on  complex  coordinate  stretching, 
which  is  shown  to  be  equivalent  to  the  anistropic  PML  medium  approach  [4]. 

In  this  paper,  we  explore  the  formulation  of  a  3-dimensional  perfectly-matched-layer  (PML)  algorithm 
that  is  used  to  describe  the  behavior  of  electromagnetic  quantities  in  outer  absorbing  boundary  layers  of  a 
nonlinear  dispersive  medium  which  serves  to  absorb  all  outgoing  waves  within  a  finite  computational 
volume.  We  consider  the  case  where  a  plane  wave  propagates  outwardly  from  a  nonlinear  dispersive 
medium  to  a  nonlinear  dispersive  PML  medium  through  a  reflectionless  PML  interface.  We  start  the 
analysis  based  on  the  extension  of  the  unsplit-field  uniaxial  PML  formulation  [4-8]  of  Maxwell's  equations 
that  are  obtained  in  the  frequency  domain  inside  the  nonlinear  dispersive  PML  medium.  We  perform  the 
inverse  Fourier  transform  of  these  equations  from  the  frequency  domain  to  the  time  domain  in  order  to 
obtain  a  set  of  ordinary  first-order  differential  equations.  Then,  these  equations  are  finite  differenced  in 
both  time  and  space  using  the  usual  staggered  Yee  FDTD  scheme  while  expanding  the  electric  and 
magnetic  field  vectors  in  time  using  the  Taylor  series  expansion  about  the  current  time  step  to  evaluate  next 
time  step  values  of  the  electromagnetic  field  quantities.  Depending  on  the  number  of  terms  kept  in  the 
Taylor  series  expansion,  we  can  numerically  update  the  field  values  to  any  desired  accuracy.  In  Section  II, 
we  use  the  piecewise-linear  approximation,  which  is  equivalent  to  using  only  the  first-order,  time- 
dependent  term  of  the  Taylor  series  expansion,  to  show  the  process  involved  in  obtaining  a  second-order 
accurate  FDTD-PML  algorithm.  To  obtain  higher-order  accurate  FDTD-PML  algorithms  in  time,  we 
simply  need  to  include  higher-order,  time-dependent  terms  in  the  Taylor  series  expansion  and  follow  the 
same  steps  shown  in  Section  II. 

In  the  absence  of  the  PML  interface,  the  FDTD-PML  algorithm  reduces  to  the  FDTD  algorithm  obtained 
for  nonlinear  dispersive  media  [9,10], 
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II.  PML  FORMULATION  FOR  NONLINEAR  DISPERSIVE  MEDIA 


For  a  wave  propagating  into  anisotropic,  uniaxial  nonlinear  dispersive  PML  media,  the  modified 
Maxwell’s  equations  under  the  PML  formulation  with  stretched  coordinates  [3]  can  be  expressed  in  the 
frequency  domain  (e10>t  convention)  as 

VxE(cd;?0  =  - ico  SPML (co) •  popu H(co; x)  ♦  (2.1) 

V  xH(co;x)  =  ico  SPML((o)#D(co;x) ,  (2.2) 

with 


Pmax 

D(a>;  x)  =  EoErE(ci);  x)  +EoV  Pp 

p=i 


(2.3) 


P J (a*  x)  =  F {  p£ (t;  x) }  =  f  {  Jdx  x£(t - x)  E(x;  x)  } =  (co)  E(co;  x) , 


(2.4) 


CO 

P^  (co;  X)  =  &  {  E(t;  x)  (t;  x)  E(t;  x)  J  dx  X^(t  -x)  [E(x;  x)  *E(x;  x)] },  (2.5) 


where  E(co;x)  is  the  electric  field  vector,  H(co;x)  is  the  magnetic  field  vector,  D(co;x)  is  the  displacement 
field  vector,  PpL(co;x)  is  the  linear  (first-order)  electric  polarization  vector,  PpNL(co;x)  is  the  nonlinear 

(third-order)  electric  polarization  vector,  SPML(co)  is  the  uniaxial  anisotropic  PML  matrix,  e0  is  the  free 

space  electric  permittivity,  eR  is  the  relative  permittivity,  jio  is  the  free-space  permeability,  pR  is  the  relative 
permeability,  XpL(t)  and  XpNL(t)  are  the  pth  terms  of  the  collection  consisting  of  time-dependent  linear 
and  nonlinear  electric  susceptibility  functions,  where  p,^  is  the  maximum  number  of  terms  which  we 
choose  to  consider  for  a  particular  formulation  of  Eq.  (2.3).  In  Eq.  (2.5)  Rp^fix)  is  introduced  to  isolate 
the  part  of  Pp^onx)  that  is  represented  by  a  convolution  function.  Also  used  in  the  above  equations  are 

notations  •  and  3F{  },  respectively,  denoting  a  dot  product  and  the  Fourier  transform  operation.  Elements  of 
the  uniaxial  anisotropic  PML  matrix,  SPML  (co),  are  given  by 


SPML(co)  = 


Sy(co)Sz(co) 

Sx(co) 

0 


V 


0 


0  0 
Sx(co)Sz(co) 

Sy(CO) 

Sx(co)Sy(co) 

Sz(co) 


(2.6) 


where  Sx(co),  Sy(co)  and  Sz(co)  are  arbitrarily  defined  co-dependent  functions  that  satisfy  the  impedance 
matching  condition  at  the  interface  of  the  non-PML  medium  and  the  PML  medium.  It  is  a  common  practice 
in  the  FDTD  community  to  choose  Sx(co),  Sy(co)  and  S2(co)  in  the  following  forms: 


* 

Sx  (co)  =  1  +  -^-  with  =  -2i_  , 

iCOEoER  £o£r  Po[Xr 


(2.7-2.8) 


2 


with 


(2.9-2.10) 


Sv(co)  =  1  +  - 


y  _ 


ICOEoEr 


£o£r  JXoJJit 


,  and 


Sz((0)  =  1+— with  — 


ICOEoEr 


EoEr  jio|iR 


(2.11-2.12) 


where  ox  ,  oy  and  az  are  the  PML  electric  conductivities,  and  at  >  ay  and  at  are  the  PML  magnetic 

conductivities  with  subscripts  x,  y  and  z  denoting  the  directions  in  which  PML  conductivities  are  assigned 
[2].  These  PML  conductivities  are  introduced  arbitrarily  in  order  to  implement  the  FDTD-PML  algorithm. 

We  first  eliminate  D(co;x)  in  favor  of  expressing  Maxwell’s  equations  in  terms  of  E(co;x),  PpL((o;x)  and 
PpNL(o);x)  by  substituting  Eq.  (2.3)  into  Eq.  (2.2).  Upon  taking  the  inverse  Fourier  transforms  of  Eqs.  (2.1), 
(2.2),  (2.4)  and  (2.5)  and  using  the  expressions  shown  in  Eqs.  (2.7)  through  (2.12),  we  obtain  the  following 
time-dependent  equations: 


no|i*To.H(t;x)+  HoM*Ti.HDelay(t;x)  +  VxE(t;x)  =  0, 
d  t  -  — 


(2.13) 


9  E(t;x)  ^  ^Ep(Lx)  d[E(t;x)Rp  (t;x)] 

EoEr . =f  +  Eo^ - ^ - +  Eo2j - 

P  P 


d  t  ^  d  t  dt 

+  eo¥o  •  [eRE(t;  x)  +  X  £p (t;  ^ + M(t;  £>  S  RpL(t;  1 

P  P 

+  eo^.[8RE^(t;x)  +  XP^lay(t;x)  +  XCMayO;2i)]-Y^H(t;x)  =  01 


(2.14) 


with 


Ep(t;x)=  JdTXp(t-x)E(x;x), 


(2.15) 


R  T  <t;x)  =  J  dx  X^(t-x)  [E(x;  x)  *E(x;  x)] , 


(2.16) 


l 

HDelay  (t;  x)  =  Jdx  £(t  -  x)  *H(x;  x) , 


(2.17) 


L 

EDelay(t;x)  =  JdxO(t-x)*E(x;x), 


(2.18) 


l 

pLDelay  =  J dT  ^  _  T)  .  pL  ^  > 


(2.19) 
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pNLDelay  (t;  g  _  J  d%  £(t  _  l)*[E(X,  x)  R  ^  (XJX)]  , 


(2.20) 


where 


(  f 


¥0  = 


,  £o£r  £o£r  £o£r 

V  / 

0 

0 


f ox  cz  °  ^ 

£o£r  £o£r  £o£r 


V 


J 


0 

f^L- 

£o£r  £o£r  £o£r 


V 


(2.21) 


°y 

\ 

°z 

£o£r 

£o£r 

J 

Eo£r 

V 

£o£r  ^ 

¥1  = 


<3>(t-x)  = 


0 

0 


exp[-  (— 2-)(t  -  x)] 

£o£r 

0 

0 


f  gx 

°y  1 

G  ^ 

Uy 

£o£r 

EoEr 

EoEr 

EoEr 

V 

J 

V 

/ 

)) 


0 

0 

£o£r  £o£r  I  £o£r  £o£r 
A 


(2.22) 


/ J 


0 

exp[-  (— ^-)(t  -  x)] 

£o£r 


exp[—  (— -z-)(t  -  x)] 

EoEr 


(2.23) 


In  the  above,  E ^(tjx),  PpLDelay(t;x)  and  PpNLDelay(t;x)  are  introduced  to  handle  the  delayed 

time-response  behavior  of  H(t;x),  E(t;x),  PpL(t;x)  and  [E(t;x)RpNL(t;x)],  respectively.  These  functions  follow 
naturally  from  taking  the  inverse  Fourier  transforms  of  convolution  functions  [l/(icoi  +  A )]  H(co;x), 

[l/(icoi  +  A)]  E(co;x),  [l/(icoi  + A )]  PpL((0;x)  and  [l/(icoi  + A )]  PpNL(co;x)  by  realizing  the  inverse  Fourier 

transform  of  [  l/(ico  1 4- A )]  is  given  by  exp(-At),  where  I  is  the  identity  matrix  and  A  is  a  time 

independent  diagonal  matrix  expressed  as  diag  [  <jx  /(e0eR),  Gy  /(Eo£r)>  gz  /(£0£r)]* 

To  solve  Eqs.  (2.13)  through  (2.20),  we  need  to  specify  expressions  for  linear  and  nonlinear  electric 
susceptibility  functions.  In  this  paper  we  consider  the  case  in  which  both  the  linear  and  nonlinear  electric 
susceptibility  functions  are  expressed  as  complex  functions  that  contain  complex  constant  coefficients  and 
exhibit  exponential  behavior  in  the  time  domain  as  follows: 


Xp  (t)  =  Re  {apexp[-(Yp)t] }  U(t) .  (2-24) 

and 

Xp1  (t)  =  Re  {ap^expRy^t] }  U(t) ,  (2.25) 
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where  Re{  }  is  used  to  represent  the  real  part  of  a  complex  function,  U(t)  is  the  unit  step  function,  and 
OCp  *  Yp  *  exp1-  and  are  complex  constant  coefficients.  Now  Eqs.  (2.15)  and  (2.16)  are  expressed  in  the 
following  forms: 


t 

Pp  (t;  x)  s  Re  {Q^  (t;x)}  =  Re  {«p  J dx  exp[-<Yp)(t-T)]  E(t;  x) } , 


(2.26) 


and 

t 

R^(t;x)  SRe{QpNL(t;x)}  =  Re{aPNL|dxexp[-(YrXt-'c)ni(x;x)*I(x;x)] }.  (2.27) 


where  complex  functions  QpL(t;x)  and  QpNL(t;x)  are  introduced  in  the  above  equations  such  that  the  real 
parts  of  these  complex  functions  result  in  PpL(t;x)  and  RpNL(t;x),  respectively. 

We  need  to  point  out  that  by  making  the  proper  choices  of  complex  constant  coefficients  and  performing 
the  Fourier  transforms  of  Eqs.  (2.24)  and  (2.25),  we  can  readily  obtain  the  familiar  constant  conductivity 

[i.e.,  otp  is  real,Yp  =  0,apL=0  and  y^=0],  Debye  [i.e.,  otp  *  Yp  .  dp1-  and  Y^  are  all  real]  and  Lorentz 

[i.e.,  ocp  and  cep1-  are  both  imaginary,  and  yp  and  y^  are  both  real]  forms  of  the  complex  permittivity  in  the 
frequency  domain. 

To  derive  FDTD  expressions  based  on  Yee  FDTD  scheme,  Eqs.  (2.13),  (2.14),  (2.25),  (2,26)  and  (2.17) 
through  (2.20)  have  to  be  solved  numerically  for  H(t;x),  E(t;x),  QpL(t;x),  Qp^tjx),  HDelay(t;x),  EDelay(t;x), 
PpLDeiay(t;x)  and  PpNLDelay(t;x)  at  each  time  step  by  correctly  carrying  out  the  numerical  integration  of 
convolution  integrals  QpL(t;x),  Qp^fix),  HDeIay(t;x),  EDelay(t;x),  PpLDelay(t;x)  and  PpNLDeIay(t;x).  Therefore, 
the  whole  solution  rests  on  the  question  of  how  to  carry  out  the  numerical  integration  of  QpL(t;x),  QpNL(t;x), 
HDelay(t;x),  EDelay(t;x),  and  PpNLDelay(t;x)  at  each  time  step.  For  that  reason,  the  rest  of  this  section 

is  devoted  to  the  numerical  formulation  that  treats  QpL(t;x),  Qp^Xtjx),  HDeIay(t;x),  EDelay(t;x),  PpLDelay(t;x) 
and  PpNLDeIay(t;x)  into  the  overall  FDTD  scheme  based  on  the  recursive  convolution  approach. 

We  first  convert  the  convolution  integrals  QpL(t;x),  QpNL(t;x),  HDelay(t;x),  EDelay(t;x),  PpLDelay(t;x)  and 
PpNLDeIay(t;x)  into  the  following  equivalent  first-order  differential  equations: 


3Qp(t;x) 

at 


+  (Yp)Qp(tii)  =  4E(t;x), 


(2.28) 


3Q^(t;x) 

at 


+  (Y^)  (t;  x)  -  apIL[E(t;x)«E(t;  x)] , 


(2.29) 


aH^djx) 

at 


+£(t).HDelay(t;x)  =  H(t;x), 


(2.30) 


aEPelay(t;x) 

at 


+  0(t)*EDelay(t;x)  =  E(t;x), 


(2.31) 


aQLDelay(t;x) 

- +£(t)*Q^elay(t;x)  =  Qp(t;x), 


(2.32) 
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+£(t).Qp—  (t;  x)  =  E(t;  x)  (t;  x) , 


(2.33) 


dQ^Delay(t;x) 

at 


where  complex  functions  0PLDelay  (t;x)  and  QpNLDday  (t;x)  are  introduced  in  Eqs.  (2.32)  and  (2.33)  such  that 
the  real  parts  of  these  complex  functions  result  in  PpLDelay(t;x)  and  PpNLDelay(t;x),  respectively. 

To  show  how  we  can  use  Eqs.  (2.13),  (2.14)  and  (2.28)  through  (2.33)  to  come  up  with  a  3-D 
FDTD-PML  algorithm  for  nonlinear  dispersive  PML  media,  we  integrate  Eqs.  (2.13)  and  (2.30)  from 
t=(n-V£)At  to  t=(n+^)At,  and  Eqs.  (2.14),  (2.28),  (2.29)  and  (2.31)  through  (2.33)  from  t=nAt  to  t=(n+l)At. 
Then  Eqs.  (2.28)  through  (2.33)  are  solved  exactly  using  the  integrating  factor  technique  for  a  given 
discrete  time  interval  to  go  forward  in  time  by  At.  The  result  is  that  we  need  to  perform  definite  integrals 
that  appear  in  the  following  equations: 


(n+‘/i)At 


(n+te)At 


P4U  f  dx  *)  +  fdxH(x;x) 

J  dx  =  J 


(n-'A)At 


(n-Vi)At 

(n+*/2)At 


(n+1/2)At 


+  JdxHDelay(x;x)  +  JdxVxE(x;x)  =  0, 


(n-‘/2)At 


(n-Vi)At 


(2.34) 


(n+l)At 


EoEr 


J  dt 


nAt 


9E(x;x) 

3x 

(n+l)At 


(n+l)At  aQL(T;x) 


(n+l)At 


+eoRe{^  J  dx  — }  +  e0Re{^  J  dt 


a  [E(x;x)  Q^(t;x)] 


P  nAt 


3x 

(n+I)At 


P  nAt 


3x 


+  EoEr^Fo  •  J  dx  E(x;  x)  +  Eo^Po  •  Re  { J  dx  Q  ^  (x;  x)  } 

nAt  P  nAt 

(n+l)At  (n+l)At 

+  eoTo.Re{^  JdTE(x;x)Q^(x;x)}  +  eoeR^.  JdxEDelay(x;x) 

P  nAt  nAt 

(n+l)At  (n+l)At 

+  eo*Fi.Re{£  JdxQj;Delay(x;x)  }+e»^i.Re{^  JdxQ^iLDelay(x;x) } 
P  nAt 

(n+l)At 

-  Jdx  V  xH(x;  x)  =  0 , 


P  nAt 


(2.35) 


nAt 


(n+l)At 

Qp  (nAt  +  At;  x)  =  exp[-(vp)At]  [  (nAt;  x)  +  otp  Jdx  exp[-(Yp)(nAt  -  x)]  E(x;  x)  ] ,  (2.36) 

nAt 


(nAt  +  At;  x)  =exp[-(y^L)At]  [  (nAt;  x) 

(n+l)At 

+  otp1"  J  dx  exp[-(Y^IL)(nAt  -  x)]  [E(x;  x)  •  E(x;  x)]  ] ,  (2.37) 

nAt 

HDelay(nAt  +  V4At;  x)  =  exp(-£At)  •  [  H1^  (nAt  -  *4At;x) 

(n+Vi)At 

+  J  dx  exp[-£(nAt  -!4At  -x)]  •  H(x;  x)  ] ,  (2.38) 

(n-*A)At 
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(n+l)At 

EDelay(nAt  +  At;x)=exp(-£At)*[EDelay(nAt;x)+  Jdx exp[-J>(nAt-x)]*E(x;x)  ] ,  (2.39) 

nAt 

(n+l)At 

q  LDeiay  +  At;x)  =  exp(-d>At)  •  [  Q  LDe'ay  (nAt;  x)  +  J  dt  exp[-0(nAt  -  x)]  •  Q  ^  (x;  x)  ] 

nAt 

(n+l)At 

=  exp(-£At)  •  [  Q^Delay  (nAt;  x)  +  Jdx  exp[-d>(nAt  -  x)  - 1  (Yp  )(t  -  nAt)]  •  (nAt;  x) 

nAt 

(n+l)At  x 

+  ap  Jdt  J dx’exp[-d>(nAt  -  t)  - 1  (Yp)(x - 1’)]  •  E(x’;  x)  ] ,  (2.40) 

nAt  nAt 


(n+l)At 

QNLE)eiay(nAt  +  At;i)  _  exp(-£At)  •  [  Q^0^’  (nAt;  x)  +  J  dx  exp[-£(nAt-x)]*E(x;x)  Q^(x;x)  ] 

nAt 

=  exp(-£At)  •  [  Q^,LDelay  (nAt;  x) 

(n+l)At 

+  J  dx  exp[-£(nAt  -  x)  - 1  (y^)(x  -  nAt)]  •  E(x;  x)  QpL  (nAt;  x) 

nAt 

(n+l)At  x 

+  o£L  J  dx  J  dx’  exp[-£(nAt-x)-I(Y^L)(x-x’)]*E(x;x)[E(x’;x).E(x’;x)]],  (2.41) 

nAt  nAt 

Furthermore,  some  of  the  definite  integrals  that  appear  in  Eqs.  (2.34)  and  (2.35)  are  manipulated  and  cast  in 
the  following  forms: 

(n+'/2)At  (n+Vi)At 

J  dx  H  Delay  (x;  x)  =  J  dx  exp[-£  (x  -  nAt  +  ViAt)]  •  HDelay  (nAt  -  V4  At;  x) 

(n-^)At  (n-'^)At 

(n+Vfc)At  x 

+  Jdx  J  dx’  exp[-4>(x  -  x’)]  •  H(x’;  x) ,  (2.42) 

(n-'/4)At  (n-H)Al 

(n+l)At  (n+l)At 

J  dx  Qp  (x;  x)  =  J  dxexp[-(Yp)(x  -  nAt)]  (nAt;  x) 

nAt  nAt 

(n+l)At  x 

+  4  Jdx  J dx’exp[-(Yp)(x-x’)]  E(x’;  x) ,  (2.43) 

nAt  nAt 

(n+l)At  (n+l)At 

J  dxE(x;x)  QpL(x;x)  =  J  dxE(x;x)  expHy^Xx-nAt)]  Q^(nAt;x) 

nAt  nAt 

(n+l)At  x 

+  apL  Jdx  J dx’exp[-(Y^L)(x-x’)]  E(x;x)  [E(x’;x)*E(x’;x)] ,  (2.44) 

nAt  nAt 
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(n+l)At  (n+l)At 

JdtEDelay(T;x)=  J  dx  exp[-d>  (x  -  nAt)]  •  E1*1^  (nAt;  x) 

nAt  nAt 

(n+l)At  x 

+  Jdx  j,dx,exp[-0(T“T,)]#E('t,;x)’  (2.45) 

nAt  nAt 

(n+l)At  (n+l)At 

J  dTQ^'^CTix)  =  J dxexp[-3>  (x  -  nAt)]  •  Q^Delay  (nAt;  x) 

nAt  nAt 

(n+l)At  x 

+  Jdx  Jdx’exp[-£x-I(Yp)(T-nAt)]«Qp(nAt;x) 

nAt  nAt 
(n+l)At  x  t’ 

+  ctp  Jdx  Jdx,Jdx"exp[-£x-I(Yp)(T-x’’)]*E(x";x),  (2.46) 

nAt  nAt  nAt 

(n+l)At  (n+l)At 

J dx  Q^LDelay  (x;  x)=  J dx  exp[-£  (x -  nAt)]  •  Q^LDelay  (nAt;  x) 

nAt  nAt 

(n+l)At  x 

+  Jdx  J dx’exp[-£  x- 1  (y^)(x-nAt)]  •  E(x’;x)  (nAt; x) 

nAt  nAt 

(n+l)At  x  x’ 

+  a£L  Jdx  J  dx’  J  dx"  exp[-<t>  x  -  I  (Y^)(x -x" )]  •  E(x’;  x)  [E(x” ;  x)  •  E(x" ;  x)] .  (2.47) 

nAt  nAt  nAt 

To  obtain  second-order  accuracy  in  time  from  a  finite  differencing  technique,  H(t;x)  and  E(t;x)  are  taken 
to  be  piecewise- linear  continuous  functions  over  the  entire  temporal  integration  range  such  that  H(t;x)  and 
E(t;x)  change  linearly  with  respect  to  time  over  given  discrete  time  step  intervals.  It  is  equivalent  to  saying 
that  we  use  only  the  first-order,  time-dependent  term  of  the  Taylor  series  expansion  for  H(t;x)  and  E(t;x), 
respectively,  that  are  expanded  in  time  about  the  current  time  step  of  t=(n-1/2)At  for  H(t;x)  and  the  current 
time  step  of  t=nAt  for  E(t;x).  Mathematically,  we  can  express  H(t;x)  and  E(t;x)  in  the  following  forms  in 

terms  of  (H)ijkn’*,  (H)ijktt+,/2,  (E)ijkn  and  (E)ijkn+1  where  superscripts  n-*4,  n,  n+'A  and  n+1  are  used  to  denote 
discrete  time  steps  at  t=(n-Vi)At,  t=nAt,  t=(n+Vi)At  and  t=(n+l)At,  respectively.  Subscripts  are  used  to 
denote  discrete  spatial  locations,  x=[iAx,  jAy,  kAz]  for  E(t;x)  and  x=[(i-^)Ax,  (j-Vi)Ay,  (k-Vi)Az]  for  H(t;x) 
with  Ax,  Ay  and  Az  being  the  spatial  grid  sizes  in  the  x,  y  and  z  directions,  respectively. 


H(t;x) 


(H)^/2  +  A  [t  -  (n  -Vi)  At]  +  higher  order  terms, 


At 


0, 


for  0  <  (n  -  Vi) At  <  t  <  (n  +  Vi) At  (2.48) 


for  t  <  0 


E(t;x)  = 


[(E)  Sv1  ~(E)Hk] 

(E),“k  + - —a - - —  (t  -  nAt)  +  higher  order  terms. 


At 


for  0  <  nAt  <  t  <  (n  + 1)  At 


(2.49) 


0, 


for  t  <  0 
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Although  we  are  not  going  to  investigate  higher  than  second-order  accuracy  in  time  in  this  paper,  it  is 
possible  to  obtain  higher-order  accurate  FDTD  algorithms  by  simply  including  more  terms  beyond  the 
first-order,  time-dependent  term  in  the  above  Taylor  series  expansion. 

Substituting  Eqs.  (2.48)  and  (2.49)  into  Eqs.  (2.34)  through  (2.47)  and  performing  the  time  integration 
from  t=(n-‘/2)At  to  t=(n+Vi)At  for  field  values  that  depend  on  the  magnetic  field  [i.e.,  H(t;x)  and  HDday(t;x)], 
and  from  t=nAt  to  t=(n+l)At  for  field  values  that  depend  on  the  electric  field  [i.e.,  E(t;x),  QpL(t;x), 
Qp^ftjx),  E^'^ftjx),  QpLDelay(t;x)  and  QpNLDelay(t;x)],  Eqs.  (2.34)  through  (2.47)  are  cast  into  the  following 

cubic  algebraic  expressions.  These  expressions  are  used  to  update  field  values  (H)pn+!*,  (E)pn+I,  (QpL)jjkD+l, 
(0,V.  (EDelay)ijkn+1,  (HDeIay)ijk"+>4,  (QpLDelay).jkn+l  and  (QpNLDelay)ijkn+1  at  each  time  step: 

Co  •QDjjk’4  +£,*(H)^  +£2  •(HDelay)uk^  +Se  =0,  (2.50) 

(Ao)i"k  *©3k+1  +(A,)i“k  •  <©*  +Re{£  r^o*(Qp)p  }+£2*(IDelay)Sk 

p 

+  Q3.Re{^(QpU5elay)1sk  } 

P  P 

+Re{^r^-}.(E)?k  [(i)?k  .(I)Kk]  +Re{^r^}*(E),"k  [(E)Sk  «(E)1“kI] 
p  -  p 

+  Re (K).(E)"k  [ ©p1  .©sk+1]  +Re  {£Cl*(E)”k+1  [ ©*  •<©*! 

P  -  P 

+  Re  [  ®)«k  •©*]  +Re  I  ©P  *©p] 

p  -  p  - 


+  Sh=0,  (2.51) 

(Qp)p1  =  ©p.0  [  (Qp)ijk  +  ©P,1  ©Sk  +©pL2  ©«k+I] .  (2-52) 

(Qp^  ),”k+1  =  ©p,o  [  (Qp^  ),"k  +  ©p^  [  ©,“k  •  ©Uk  ] + ©P^  [  ©,“k  •  ©3k 1  ] 

+  0p^[©3k+1*©Sk+I]].  (2-53) 

(HDe,ay)Sr  =  £4  *[  (HDeIay)Sk^  +£5  •<H)j*4  +£5  •0Djjk%]  •  (2-54) 

(EDelay)p1  =£4 ' •[  dDehy)|Jk  +£5  *©3k  +£6*©3k1],  (2.55) 

(Qray)p  =£4*[(QjDelay)3k  +np.o*(Qp)3k  +n^,*(E)3k  +n^.©Hk+1],  (2.56) 


(QpNLDe,ay)3k+1  =£4  •[(QpNLDelay)3k  +(Q?%  nM  *©3k  +(QpNL)3k  Hpa-  *©3k’ 

+  •  (E)p  [  (E)Jk .  (E)p] + n£  •  ©*  [  ©nk .  (e)^+i] 

+n^*(E),nk  [  <©£*  •(©;k+1]+ng.(E)nk+1  [  (E)p  *  ©3k  ] 

+£$  •(E)j’k+1  [  (E)3k  *(E)|Sk+,]  +  n^.(E)Hk+1  [  ©3k+1*(E)^]],  (2.57) 

where  SE  and  SH  are  given  by 
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^^[(EZ),W  - (E z )"(H«k  ] -  (JioMR)Az  t(E y >>V+H)  -(Ey).J(k^)] 

JU  [(E  X  )^+V4)  " (Ex  )*k"^>  ]  ~  Ohu)Ax  *-(Ez  )”,+,/4)jk  _(Ez)”i_V4)jk-1 

J)A^yW  “ (E y >G-WP ] - (MiR)Ay  t(E x ).W  -(Ex).Wl  I 


'(eoeR)Ay[(Hz)^,k  (Hz  ]^>k  ^ +  (e„br)Az  ^(H  y  (Hy)i^^ 


and  (A0)!k  and  (A,),“k  are  the  matrices  that  depend  on  (Q^lp  ,  material  properties  and  At,  and  0pO, 


0L 

UP,1’ 

0L 
UP,2  * 

Up,0  * 

oNL 
Up,l  * 

Up,2  * 

e1^ 

WP,3  » 

Qq  » 

||P 

VP 

,  Q3 , 

Q4  9 

Q5 , 

n  -  ■ 

rNL 
’  1  p,0 

pNL  pNL 

’  A  p,l  9  1  p,2  * 

l-NL 
AP,3  i 

rNL 
►  Ap,4  - 

T->NL 
*  Ap,5  ■ 

nL 

>  AAp,0 

nL 

’  AAp,l 

’  AAp,2 

TT  NL 
>AAp,0 

»  A1p,i  » 

FfNL 
AAp,2  1 

T-rNL 

»  1Ap,3 

j  AAp,4  9 

nNL 

1  Ap,5  * 

AAp,6 

and  fl$  are 

the  coefficients  and  matrices  that  depend  only  on  material  properties  and  At.  The  material  properties 
required  for  the  evaluation  of  these  coefficients  and  matrices  are  ap ,  Yp  »  otp’L  >  YpL »  ox  >  oy  and  oz  ■ 

Shown  in  Appendix  are  the  expressions  of  these  coefficients  and  matrices. 

Using  the  above  FDTD-PML  algorithm  the  computer  simulation  can  be  performed  for  electromagnetic 
waves  that  propagate  inside  nonlinear  dispersive  PML  media  by  simply  going  through  the  following  steps: 

(1)  First,  as  part  of  the  initial  condition,  time-invariant  coefficients  0p,o*  ®p,i>  ®p,2  *  ®p!o  *  * 

and  0^3 ,  and  time-invariant  matrices  Q0  ,  ,  Q2  ,  Q3  ,  Q4  ,  Q5  *  Q6  ,  r£0  *  , 


•pNL  rNL  rNL  rNL  pNL  i-rL  ttL  t-tL  ttNL  ttNL  ttNL  ttNL,  ttNL  ttNL 
f  p,l  *  fp,2  *  A  p,3  »  A  p,4  »  A  p,5  *  AAp,0  »  AAp,l  9  AAp,2  ’AAp,0  9  AAp,l  9  AAp,2  *  AAp,3  *  AAp,4  >  lxp,5  > 

and  11$  are  all  calculated  at  the  beginning  of  the  simulation  for  given  values  of  otp  .  Yp  > 


ap1- *  YpL *  ap  .  ax  >  oy ,  az  and  At.  These  values  are  stored  in  computer  memory  and  used  in 
calculating  the  updated  field  values  at  each  time  step. 

(2)  Using  Eq.  (2.50),  (H)ijkn+H  is  calculated  based  on  the  known  values  of  (Hip"'**  and  (HDelay)ijkn'v'2  and 
(E)i/. 

(3)  Using  Eq.  (2.54),  (HDelay)ijkn+/;  is  calculated  based  on  the  known  values  of  (HDelay)ijkn'!/2,  (H)ijkIMr2 
and  (H)ijk“  ^  • 

(4)  Using  Eqs.  (A.8)  and  (A.  12),  (A0),"k  and  (Aj  ),“k  are  calculated,  respectively,  with  the  values  of 

(Q»“V- 

(5)  Using  Eq.  (2.51),  (E)ijk“+1  is  calculated  based  on  the  known  values  of  (E)ijkn,  (E  y)ijkn,  (£>p  )ykn, 
(QpLDelay)ijkn  and  (H)ijkn+/2.  Because  Eq.  (2.51)  represents  coupled,  cubic,  algebraic  equations,  the 
nonlinear  Newton-Raphson  method  [11]  is  used  to  solve  for  (E)ijkn+I  by  finding  the  zeroes  of 
Eq.  (2.51)  with  (E)jjkn  as  the  initial  guess. 

(6)  Using  Eqs.  (2.52),  (2.53),  (2.55),  (2.56)  and  (2.57),  (QPL) ijkn+1,  (Qp^V*1,  (EDelay)ijkn+l, 
(OpLDday)ijtn*1  and  (QpNLDelay)ijkn+1  are  calculated  based  on  the  known  values  of  (E)ijkn+I,  (E)ijk“, 

0QpL)ijk”,  (Qp“0«k".  (B6*9)*".  (QPLDe,ay)ijkn  and  (QpNLDelay).jkn . 

(7)  Increment  the  time  step  by  At.  Go  back  to  step  (2)  and  repeat  the  whole  process  over  again. 

Shown  in  Figure  1  is  the  flow  chart  of  numerical  steps  required  to  update  field  values  as  described  above. 
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In  the  absence  of  the  PML  interface,  ax  *  cry  and  oz  are  all  set  to  zeroes  and  the  nonlinear  dispersive 

FDTD-PML  algorithm  reduces  to  the  usual  nonlinear  dispersive  FDTD  algorithm  obtained  for  a  nonlinear 
dispersive  medium  with  the  third-order  electric  susceptibility  functions  [9,10]. 

Also,  we  get  the  linear  FDTD-PML  algorithm  by  simply  setting  otp1'  to  zero.  Furthermore,  if  both  ap 

and  a  p1^  are  set  to  zeroes,  the  FDTD-PML  algorithm  reduces  to  the  case  of  the  simple  PML  algorithm  in 
free  space. 

III.  CONCLUSIONS 

We  present  in  this  paper  the  formulation  of  a  three-dimensional  FDTD-PML  algorithm  inside  nonlinear 
dispersive  PML  media  that  is  used  to  absorb  all  outgoing  electromagnetic  waves  within  a  finite  simulation 
volume  to  create  the  notion  of  infinity  at  the  outer  layer  boundary  of  the  computational  volume.  Because  of 
the  use  of  the  piecewise-linear  approximation,  the  FDTD-PML  algorithm  provides  second-order  accuracy 
in  time  for  the  calculation  of  electromagnetic  field  quantities.  The  resulting  forms  of  the  FDTD-PML 
algorithm  tell  us  that  we  need  to  solve  coupled  cubic  equations  for  the  three  components  of  the  electric 
field  vector  at  each  time  step  due  to  the  nonlinear  behavior  of  the  third-order  electric  susceptibility 
function.  In  the  absence  of  the  nonlinear  contribution,  the  electric  field  components  are  no  longer  coupled 
and  each  component  of  the  electric  field  vector  can  be  updated  independently  of  the  other  two  components. 
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APPENDIX 


This  appendix  gives  the  explicit  expressions  of  coefficients  and  matrices  seen  in  Eqs.  (2.50)  through 
(2.57).  The  coefficients  and  matrices  shown  below  are  for  ©p0,  ©p,i,  ©p,2»  ©p!o>  -  ®pl2-  ®p!s  > 

(Aq  ) jji;  ,  (A)  ,  O0  ,  ,  02  ,  O3  ,  O 4 , 05 ,  £2g  ,  f  po  ,  Tp^  ,  Tpj  ,  Fp  2  >  rp3  ,  rp4  ,  Fp  5  , 

Q  ,  n£j,  n^2  ,n$  ,  nJJ  ,  n$  ,  n£  ,  n£  .  n£  ,  and  ng.  Also,  to  express  these 

coefficients  and  matrices  in  more  compact  forms,  additional  terms,  such  as  qp  0 ,  qp vl ,  ^p,2  *  Sp,o  *  Sp,i  * 
O*  Sp.l  -  Sp,2-  9pj’  <0  ?M’  (Vo)*>  (Vl)x»  (V2)x>  (<Pl)x.  (<p2>x»  (Cp.oL-  (Cp,l)x  > 

Ox  .  (Ox  .  (Ox  -  (Ox  -  (Ox  -  (Ox  -  (Ox  .  (Ox  -  Ox  •  Ox  -  Ox  -  Ox  . 

(Ox  -  Ox-  Ox-  (Ox-  (Ox-  (Ox-  (Oxand  (Ox*  are  defined-  These  additional 

terms  are  shown  following  the  expressions  for  coefficients  and  matrices. 


0p,o=O 

0p,i  =^rKk-0. 


£r 


qI  -  aP  gL 

Op,2  -  £r  4p,2- 

©NL  _  fc  NL 
p,0  ~  Sp,0  * 


\ NL  _  aP  rftNL 


©NL  _  r&  NL  o  £  NL  ,  fc  NL 

p,l  “ - LSp,  1  ~ZSp,2  +Sp,3i’ 

£r 


©p^=^[20-0- 

£r 


„NL 

oNL  _  aP  £  NL 

Op,3-  — ^p,3. 


(A0),"k  = 


[<A0>Sklii  0 

0  [(A0)|jic  ]22 


0 

0 


0 


0  [(A0)Sk]33 


(A.1) 

(A.2) 

(A.3) 

(A.4) 

(A.5) 

(A.6) 

(A.7) 

with  three  diagonal  elements  expressed  as  (A.8) 
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[<Ao  ),"k  ]„  =  1 +[(—)  +  O)  -  (— )]  v + [(— )  -  (-^)][(— )  -  (— )1  (cp2)  > 

£o8r  £o£r  £o£r  2  £o£r  £o£r  £o£r  £o£r 


+— Re{ JaHp,o^P.i  }+— t(7f)  +  (7r)-(7r)]Re{SaP^p.2  j 

Er  •*“ '  Er  EoEr  EoEr  EoEr  „ 

p  p 

+ — [(— ^-) - (— ^— )] [ (— ^-)  - (-*-)] Re  O*  } 

£r  8o£r  £o£r  EoEr  £o£r  p 

+fRe(S?S<«pNI'>5«  l 

Er  p 

+-[(—)  +  (—)-(— )]Re{V^  (<2?%  } 

8r  £o£r  £o£r  £o£r  "  k  J 

+-!-[A)-(-^)][<i-)-(-^)]Re  {Y  <£((5?% ) , 

Pd  PftPn  PdPd  PoPd  P«PD  r  J 


KAo).jk]22  =  Replace  [x->y,  y->z,  and  z->x)  in  [(A0)ijk Li  > 
[(A0)yk  ]33  =  Replace  [x-»z,  y->x,  and  z-»y]  in  [(A0),“k  ],, , 


(A.9) 

(A.  10) 
(A.  11) 


(Ai)Sk  = 


[(A,)J^]„  0  0 

0  [(A,  >3^22  o 

o  o  [(A,  )yk  ]33 


with  three  diagonal  elements  expressed  as  (A.  12) 


[(AOjUj  =  - 1  +  [  A~)  +  O)  - O)] ^  +  [(— )  - (—)][(—)  -  (— )] [((P, ) x  -(cp2)xl 


EoEr  EoEr  EoEr  2  EoEr  EoEr  EoEr  EoEr 

1  Ov 

— [(— . 

£r  £o£r  £o£r  £o£r 


P  P 

+-[(^)-0)][0)-(-^)]Re{Y  c*W£,>*  -0x1} 

£r  Eo£r  £o£r  £o£r  £o£r 

--Re{Y(Q«L)3k  }+-[0)+0)-0)]Re{X[C-^](Qr)ijk } 

Er  "  p  J  Er  EoEr  EoEr  EoEr  „ 

P  '  P 

+—[(—)  -  0)][0)  -  0)1  Re  {£  sS  (QpL  ),"* } .  (A- 13) 

Pd  PdPo  PaPd  PaPd  PdPd  K’  r  J 


[(A,  ),"k  ]22  -  Replace  [x->y,  y->z,  and  z->x]  in  [(A,  ),“k  ]u , 
[(A,  ){jk  ]33  =  Replace  [x-»z,  y->x,  and  z->y]  in  [(Aj  )j]k  ]n , 


Qn  — 


(O0)u  ® 


0 


0  (^0)22  0 

0  0  (£20)33 


with  three  diagonal  elements  expressed  as 


(A.  14) 
(A.  15) 

(A.  16) 
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(A.  17) 


(Q0)n  =l+[(— )  +  (—)"(— )]-^-+[(—)-(—)][(-^-)-(-^)](cp2)x , 

£o£r  £o£r  £o£r  2  £o£r  £o£r  £o£r  £o£r 


(Q0)22  =  Replace  [x— >y,  y-*z,  and  Z“-»x]  in  (£20 )n  » 
(£20)33  =  Replace  [x— >z,  y-*x,  and  z— >y]  in  (Q0)u , 


Qj  — 


(Qi)u  0 


0 


0  (Qj)22  0 

0  0  (^1)33 


with  three  diagonal  elements  expressed  as 


(Oi)n  =-!+[(—)+(—)-(— ~(<P2), 

£o£r  £o£r  £o£r  2  £o£r  8o£r  £o£r  £o£r 

(QO22  =  Replace  [x-»y,  y->z,  and  z^x]  in  (i^1)11 , 

(£*2)33  =Replace  [x-^z,  y->x,  and  z->y]  in  (&i)n  , 


Q2  — 


(Q2)n  0 


0 


0  (^2)22  0 

0  0  (Q  2 )  33 


with  three  diagonal  elements  expressed  as 


(Q2)ii  )](¥i),. 

£o£r  £o£r  £o£r  £o£r 


(Q2)22  =  Replace  [x-*y,  y— »z,  and  z-^x]  in  (Q2)u , 
(^2)33  =  Replace  [x— >z,  y— >x,  and  z-»y]  in  (Q2)n , 

(  1 


Q3  = 


£r 


’(^2)11  0  0 

0  — (q2)22  0 

£r 

0  0  —  (Q2)33 

£r 


with  three  diagonal  elements  expressed  as 


— 


(Vo>*  0 


0 


Qs  = 


0  (Vo>y  0 

0  0  (\|/0)z 


[(V,)x-(V2)xl  0 

0  [(Vl)y  ~(V2)y] 


0 

0 


0 


0  [(Vi)z-(V2)J 


(A.  18) 
(A.  19) 

(A.20) 

] ,  (A.21) 

(A.22) 

(A.23) 

(A.24) 

(A.25) 

(A.26) 

(A.27) 

(A.28) 

(A.29) 

(A.30) 
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(V2)x 


rL  — 

1p,0 


0  0 
0  (V2)y  0 

0  0  (Vj)z 


(rpLo),i  0  0 


(A.31) 


0  (rDLo) 


p,0  /  22 


with  three  diagonal  elements  expressed  as 


(A.32) 


0  ^  (rp,0)33 

X 

(rpL0)n  =  [  -«+[  (-)+(—)-(—)  ]  & 

r’  r’  e.o*  f'  -  r'  r- . 


£o£r  EoEr  BoBr 


+  [(— )-(—)][(—)-(—)]  (Cp,o)x  . 

BoBr  BoEr  £o£r  EoEr 


(rpLo)22  =  Replace  [x->y,  y->z,  and  z^x]  in  (rpL0)n , 
(rpLo)33  =  Replace  [x->z,  y-»x,  and  z-»y]  in  (rpL0)u , 


(A.33) 

(A.34) 

(A.35) 


nNL  _ 

Lp,0 


(  ryNL 

aP  /t>NL  \ 

U  p,o  )n 
Er 

0 
0 


0 

nNL 

aP  /t-'NL  \ 
“  V  p,0  >  22 
Er 


0 

0 

aP  /t*NL  \ 

“U  p’°  '33 
Er 


with  three  diagonal  elements  expressed  as  (A.36) 


J 


/pNLv  r  /  y  Z  \  /  X  \  ir  -NL  INL  -i 

(rp.o  )n  =  [  (— )  +  (— )  -  (— )  1  [  Sp,i  -  3 <;p,2  +  3 Sp.3  - Sp,4  ] 
EoEr  EoEr  EoEr 


+  [(!LH%iLH%Os  -2(Ox  +  Ox  "Ox -2©,+©,]. 
EoEr  EoEr  EoEr  EoEr 

(A.37) 
(A.38) 

(A.39) 


(rpNJ-)22  =  Replace  [x->y,  y->z,  and  z->x]  in  (r”o  )n  , 


NLx 


(rpN0L)33  =  Replace  [x->z,  y-»x,  and  z->y]  in  (r£o )„ 


NLx 


rNL  _ 
XP.l  - 


f  fyNL 

<*P  /T,NLx 
—  Up.1  hi 
Er 

0 
0 


0 

UP  /-p  NL  \ 
— -Up,i  )12 
Er 

0 


0 

0 

„NL 

aP 

T“^P,1  '33 
Er 


with  three  diagonal  elements  expressed  as  (A.40) 
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(rpN,L)n  =  [  A-)  +  (^)-(— )  ]  1 2$-4<$+2$  ] 

P  EoEr  £o£r  £o£r  W>  K  H 

+ [(— )  -  (—)][(—)  -  (— )]  1 2  (C^)  x -20, -2Q,  +20x3*  <A'41) 

£o£r  £o£r  £o£r  Eo£r 


(rp^)22  =  Replace  [x-»y,  y— >z,  and  z->x]  in  (r^L  )n , 
(fp^  )33  =  Replace  [x->z,  y-*x,  and  z— >y]  in  )u , 


(A.42) 

(A.43) 


■p»NL  __ 

AP,2  - 


(  ^NL 

<*P  /rNLx 
U  p,2  )ll 
£r 

o 

0 


0 

“p1”  /pNL  \ 

U  p,2  '22 
£r 


0 

0 

^NL 

£r 


with  three  diagonal  elements  expressed  as  ( A.44) 


(rpN2L)n  =  [  (— )+(^)-(— )  ]  [  ?p,3  “C  ] 

P’  £o£r  £o£r  £o£r  m 

£o£r  EoEr  £o£r  £o£r 

(T$)22  =  Replace  [x-»y,  y^z,  and  z-»x]  in  (r^  )n  * 
(rpN^)33  =  Replace  [x-»z,  y-*x,  and  z-»y]  in  (T^  )n  , 


(A.45) 

(A.46) 

(A.47) 


f  „nl 
aP  /rNLx 
V1  p,3  >>11 


r-NL  _ 
ap,3  “ 


£r 


0 

0 


0 

„NL 

aP  /t»NL  \ 

p,3  '22 
£r 

0 


0 

0 

^-<r; s-)» 

Er 


with  three  diagonal  elements  expressed  as  (A.48) 


(rp^)„=^t^-2^+^]+[(^)H^)-(^)][^-2^+^] 

£o£r  EoEr  EoEr 

+[(^L)_(^)][(^)-(J^)][0x -20.+0,], 

EoEr  EoEr  EoEr  EoEr 

(rPN3')22  =  Replace  [x->y,  y-»z,  and  z->x]  in  (r^),] , 

(rpN^)33  =  Replace  [x->z,  y->x,  and  z-4y]  in  (Fp^ )u , 


(A.49) 

(A.50) 

(A.51) 
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-p  NL  _ 

Ap,4  ~ 


(  „NL 

^-On 

£r 

0 
0 


0 

nNL 

aP  /'pNL\ 
“7  p,4  >>22 

£r 

0 


0 

0 

rvNL 

aP  /p»NL  \ 
— —  U  pf4  '33 
£r 


with  three  diagonal  elements  expressed  as  (A.52) 


£o£r  £o£r  £o£r 


+ [(— )  -  (—)][(—)  -  (—)}  [  2  ®  x  -  2  ®  x  ] 
£o£r  £o£r  £o£r  £o£r 


(A.53) 


(rpN4  )22  =  Replace  [x->y,  y->z,  and  z->x]  in  (F^  )„  , 


(rM  )33  =  Replace  [x->z,  y-»x,  and  z-»y]  in  (I^  )u 


NLn 


(A.54) 

(A.55) 


(  ~NL 

aP  /r-.NL  \ 
1^,5  211 


rNL  — 
1P-5  “ 


£r 


0 

0 


nNL 

aP  /pNL\ 
“7“^A  P>5  >22 
£r 


0 

0 

„NL 

aP  /p»NL  \ 
P.5  >>33 
£r 


with  three  diagonal  elements  expressed  as  (A.56) 


/rNL  X  _  tNL  tNL  ,  r  X  ,  ,^7.  X  /9xx  A/NLx  .NL 

(pp,5  hi  =  £p,o  ^p,3  +  [  (— )  +  (— ) - (— ) “ (Yp  )  ]  Sp,4 

£o£r  £o£r  £o£r 


+  [(— )  - (—)][(—)  - (— )]  (Ox  > 


£o£r  £o£r  £o£r  £o£r 


(A.57) 


(r^)22  =  Replace  [x— »y,  y— 3Z,  and  z-»x]  in  (r”^)u 


NL  \ 


(rp^  )33  =  Replace  [x->z,  y->x,  and  z-»y]  in  (rp"j  )u  , 


NL  i 


(A.58) 

(A.59) 


^  ap 
£r 


ttl  = 

11p,0 


(7tp,0) > 


0 

0 


ap 

£r 


(Jtp.o) ) 


0 

0 


0  t-K  »). 

£r 


(A.60) 


f  „L 

-(<2)xi 


n 


L  _ 

P.> 


£r 


0 

0 


0 

^[(<,)y-(<,2)y] 


0 

0 


•7^-[(7tp,i)z  -(Jtp,2)2] 
£r 


(A.61) 
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(  fyL 

aP 

Er 


rip, 2  = 


(O 


■p.2'x 


0 

0 


0 

ap,_L ) 

Er  (jlp-2)y 


0 

0 


0  —  (<2)z 

Er 


) 


T-jNL  = 

AAp,0 


NL 


Op_ 

Er 


[(Ox  “(Ox  I 
o 
o 


Er 


0 

l(Oy  “(Oy] 

0 


0 

0 


NL 


■^-[(Oz“(Ox] 


(  „NL 
UP  /'_NL\ 

SR  (^)X 


n1^  = 

AAp,l 


^NL 

0  ^<C>, 


o 


o 


0 

0 

„NL 

Er 


TT  NL  — 
AAp,2 


f  NL 

— (nJS)n 

Er 

o 

o 


o 

„NL 

— (O* 

Er 

o 


o 

0 

„NL 

^(n^)33 

Er 


J 


(n^)„  =(Ox  -2  (Ox  +(Ox  “(Ox  +2  (Ox  “(Ox* 

(n^2)22  =  Replace  [x-yy]  in  (II^),,  , 

(11^2)33  =  Replace  [x->z]  in  (n^)„  , 


t-tNL  — 
AAp,3 


(  nNL 
Er 

0 

0 


0 


0 


„NL 

— (np^)22  0 

Er 


0 


apL 


ffiS) 


Er 

<n$)„  =(Ox  -2 (Ox  “2(Ox  +(0 


P,3  >  33 


(np  3)22  =  Replace  [x-»y]  in  , 


(A.62) 


(A.63) 


(A.64) 


(A.65) 

(A.66) 

(A.67) 

(A.68) 

(A.69) 

(A.70) 

(A.71) 
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(A.72) 


(n$)33  =  Replace  [x->z]  in  (n£)n , 

~NL 

— (0„  0  0 

8r 

o  — (n^)22  o 

Er 

0  o  —  (n^)33 

Er 

(n^)22  =  Replace  [x-»y]  in  (n£)u  , 

(n^)33  =  Replace  [x->z]  in  (II$)n  , 

*-(n»)„  o  o 

Er 

o  — (n^)22  o 

Er 

o  o  —  (n^)33 

Er 

(npN5)u  =(<k)x  -2(^)x  +(<V)x . 

(n^5  >22  =  Replace  [x->y]  in  (n^)H  , 

(np^)33  =  Replace  [x-»z]  in  (n£)u  , 

— (n$)„  0  0 

Er 

o  —  (n^)22  o 

Er 

o  o  — (n^)33 

Er 

(n^6)ll  =  2(jl$)x  -2(7Tp,7)x, 

(n^)22  =  Replace  [x-»y]  in  (II$)n , 

(n^6>33  =  Replace  [x-»z]  in  (nj$)n  , 


(A.73) 

(A.74) 

(A.75) 

(A.76) 

(A.77) 

(A.78) 

(A.79) 

(A.80) 

(A.81) 

(A.82) 

(A.83) 

(A.84) 
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T-rNL 

AAp,7 


(  ^NL 

2=-  <*$>, 

£r 

o 

0 


0 


0 


„NL 

0 

„NL 

£r 


with  elements  defined  as  follows 
£p,osexp[-(Yp)At], 

.  r  i  l-exp[-<Yp)At] 

£p,i  =  J  dt  exp[-  (Yp)x]  =  At  [ - — - ] , 


0 

At 


(Yp)At 


^^exPHYpNL)At], 


At 

&■-  Jdxexp[-(YpNL)x]  =  At[ 


0 

At 


l-exp[-(YpL)Atl 

(YpNL)At  J> 


« ■ ■  b  £  >  - “^e*i  1~7vtff~-expHC)A1]]  • 


0 

At 


Jdx(^)2  exp[-(YpL)x] 

0 

=  -expf-CYp^At]]  -exp[-(YpNL)At][l  , 


•(YDNL)AtL[(YpNL)Atf 


ijp  jS  J*dx  Jdx’exp[-  (YP)(T  -  x’)l  =  (At) 2  — [l  - 


1  r  l-exp[-(Yp)At] 


o  o 


(Yo)At 


(Yo)At 


] 


Ar  f  x’  I  ,  1  ri  1  r  l-exp[-(Yp)At]11 

Jd^ j^’ (^r > “pmvpXt-to: = <At)  _[r_D — lJ- 


(Yp)At  2  (Yp)At 

At  t  i  r 

$ s  J ' dT  J  dT’  exp[-(Y^)(T  -  t’)]  =  (At) 2  [l  - 


0  0 

At  T 


l-expHY^At] 


0  0 


(YoL>At 


1 


(A.85) 

(A.86) 

(A.87) 

(A.88) 

(A.89) 

(A.90) 

(A.91) 

(A.92) 

(A.93) 

(A.94) 

(A.95) 
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(A.96) 


At  x 


s£  s  JdT  |  df  (■£■ )  exp[-(Y^)(T-T’)] 


0  0 


=  (At)2 


[i _ 1 _ [1- 

✓  NLn  a  ,  L 


1  ri  1  r.  l-expHy^At] 


(Yo^At  *"2  (YoNL)At 


(Y0NL)At 


■]]• 


At  t 


^  =  J ' dT  J  dT’  ) 2  exP[~(Y^L)(T  -  x’)l 


0  0 

•  2  i  r  i  i  2 

=  (At)2 — xsr — I  - - +- 


(YP  )At‘-3  (y"L)At  [(Y"L)Atf  [(Yp  )At] 


-[l-expHyl^At]]], 


At  x 


Jdx  J  dx’  ) 3  exp[-(Yp  L)(x  -  x’)] 


0  0 

2 


-  (At) 


-}—[i _ 1 

1,NL\A*1-  a  /..nl 


6  6 
+  - 


(Yp)AtL  4  (Yp  )At  [(Yp  )At]2  [(Y^)At]3  [(Yp  )At] 

For  the  x  component: 

(V0)x  =  exp[-(-^-)At] , 

£o8r 


x[l  -  exp[-(Y^L)At]]] » 


At  1  — expH— — )At] 

(Vi)x  h  f dx exp[-( — — )x]  =  At [ - ^ - J, 

J  £o£r  .  aY 


( — ~)At 

£o£r 


^  1-expK— ^)At)] 

(v2)  x  =  f  dT  (~ )  exp[-(— )x]  =  At - [ - — - exp[— ( — — )At)]] , 

J  At  £o£r  ,  ,  A  .  ,  O’y  v  x  £o£r 


(— )At  ( — — )At 

£o£r  £o£r 


At  T 


-  ■  ,  r  l-exp[-(— 2-)At]1 

(<P,)x  -  [<h  fdT’exp[-( — — )(t-T’)]  =  (At)2 - Ll - — - J, 

£o£r  £o£r 


0  0 


At  T 


«p2)x  =  fdx  fdx’(^)expH— Xx-t’)] 

J  J  At  EoEr 


0  0 


,  ri  ,  1-expK— 5-)At] 

: (At)2  - —  [  1 - ^ - ]J  , 

Ox  L2  Cfv  ax  j 

(— -)At  (— )At  (— )At 

EoEr  £o£r  £o£r 


(A.97) 


(A.98) 


(A.99) 


(A.  100) 


(A.  101) 


(A.  102) 


(A.  103) 
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(A.  104) 


At  T 


LXl  t 

(Cp,0)x  s  JdT  J dx’exp[-(Yp)T’]  exP[_(“)(T_T’)^ 


0  0 


=  (At)2 


-[ 


l-exp[-(Yp)At] 


£o8r 

At  T  T’ 


[(— )-(Yp)]At  (YP)At  (— )At 

K  £o8r 


(Cp,,)x  =  Jdt  J dT’J dT"  exP[~(Yp)(T’_T")]  expH-^-XT-'c’)] 


0  0  0 


=  (At)3  — J“ — {  [~ ~ — ][l- 

(Yp)At*  «- 


l-exp[-( — — )At] 


£o£r 


(~2_)At 

£o£r 


(— )At 

£o£r 


] 


-][ 


l-«pHYp)Al]  1_eKph(li,4tl 


( — ^)At-(Yp  )At 

£o£r 


(Yp)At 


ov 

(— )At 

£o£r 


]}• 


At  T  T’  H  Q 

(Cp,2)  x  =  JdT  J  dT’ J  dx"  )  exp[-(Yp)(T’-T" )]  exp[-(^-)(x  -  t’)] 


o  o 


i  rr  1  -11-1  l-exp[-(— *-)At]  , 

-(A03  1  {  [ _ l _ 1  fi — fi _ §!** - ][ - 1 - +—!—]] 

-(At)  „  JL9  LI  CT  Jl(o^)At  (Y^)AtJJ 

£o£r 


(Y^At1^,*,  *  (i_)4t 

£o£r  £o£r 


4 


1 


ir  i  ~i r*~exp[~(Yp)At]  1-CTpHi±)A,1-|  \ 


(i)4,.(Yx)4t  (Yp)AI  <Y> 

£o£r 


(-^)At 

£o£r 


At  T 


(Qx  *  JdT  J dT’expt-CY^)!’]  exp[-(-^-)(T-T’)] 


0  0 


=  (At) 2 


rl-expKY.^At] 

L  /..NLxaa  rr  J’ 


[(-^)a.-(y5"W  <Y^)At 

EoEr 


(— *~)At 

EoEr 


(A.  105) 


(A.  106) 


(A.  107) 
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At  x 


0*s  fdx  f dx’(f ) exp[-(YpL)x’] exp[-(-^4(x-x’)] 

J  J  At  K  Eo£r 


0  0 


'■  (At) 


•{[l-Bp(-«>]W 


(^)At-(Y"L)At 
EoEr  k 

-[l-exp[-(Y^)At]][ 


[(Y^  )At]2 


1 


EoEr 


_ex±^fl+[i-e  „  )4t]][ 


1 


OAt 


£o£r 


[(-^)At][(-^)At-(Y^)At] 

£o£r  £o£r 


]}. 


At  t 


(©,  =  Jdx  j* dx’J* dt"  exp[-(YpL)(x-x")]  expH^M(x-x’)] 


0  0 


1  rr  1  ir  l-exp[-(— )At]_. 

-  <">Jrsnr{  [-rM  [i - =->»—] 


(Yp^)At  (^)At- 
£o£r 

i 


(— )At 

£o£r 


'(~^~)At  -  (Yp^  )At 

£o£r  h 


][ 


l-exp[-(Y?L)4t)  ‘-“PHrr)41! 


tCw 


_ £o£r 

av 

(— )At 

£o£r 


]}. 


At  X  X’ 


Ox  s  fdx  f dt’  f dt" )  exp[-(Yp)(x’-x")]  exp[-(— )(x-x’)] 

K’  J  J  J  At  F  £o£r 


0  0 


1  rr  1  -in  l-exp[-(— *-)At] 


— — ]] 

NLxa.JJ 


(YPNL)AtlL(^)AtJL2  (^)At  (^)At  (YpL)At 


£o£r 


£o£r 


EoEr 


[ 


1  irl-expKYp^M] 


]  [ — 1 — ]  [ 

JL/.,NLaa*JL 


(-^-)At-(Y^IL)At  ^ p  ^At 


EoEr 


(YpWL)At  (^)At 

EoEr 


•j: 


(A.  108) 


(A.  109) 


,  (A.  110) 
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1 


At  T  T’  n 

O*  s  JdT  J dt’J dT"  )2  exP[-(YpL)(t-t")]  exp[-(^)(T-T’)] 


0  0  0 


=  (At)3 


-{±-— ‘ 

^  Q  rr 


[(Yp1*  )At][(— i-)At]  3  (— )At  [(— ^-)At]2  (YpL)At 


EoEr 


EoEr 


EoEr 


[(Yf  )At][(— )At] 


EoEr 


[l-exp[-(-^-)At]][ 


EoEr 


[(— ^~)At]3  [(Y^)At][(^)At]J 

£o£r  £o£r 


[(Y^)At]2[(^)At]  [(YpL)At][(— )At  -(Y^)At] 


-[l  -exp[-(Y^IL  )At]]  [- 


8o£r 

2  [(— -5— )At] 

£o£r 


£o£r 


[(Yp1')At]3[(^)At-(Y^L)At] 


-]}. 


(C^)x  s  JdT  } dx’ j dT"  )  exp[-(YpL)(t-x")]  expH%-T’)] 

0  0  0 

=  (At)3 


[(Y  p  L )  Atl  [( — ~)  At  ] 

£o£r  Bo£r 


{- - - — 

2  (~i_)At 


+  [l  -  exp[-(-^-)  At]  ]  [• 


EoEr 


[(—*-)  At] 2  [( — — )At  -  (Y^1-  )At]2 


£o£r 


-[l-exp[-(YpL)At]]  [- 


£o£r 

(— *-)At 

£o£r 


[(Yr)At]"[(— )At-(Yr>At] 


NL> 


£o£r 


(— *-)At 

EoEr 


[(Y  p1' ) At]  [(— i-) At  -  (y” L ) At] 2 

EoEr 


-exp[-(Y^)At][- 


(— )At 

EoEr 


[(Y?L)At][(— )At-(Yp")At] 


NL  \ 


-]}. 


EoEr 


(A.  Ill) 


(A.  112) 
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At  T  T 


®x  =  Jdx  J  dx’j  dx"  (■£  )(-£ )  exp[-(YpL)(t-x" )]  exp[-(-^-)(x-x’)] 


0  0  0 


=  (At)! 


[(Y^)At][(-^-)At]  3  (-^)At  [(— )At]2  [(Y^)At][(— )At] 

8o£r  £oEr  £o£r  £o£r 

-[l_exp[-(^-)At]]  [— ^ - + - - - 


t(— ^~)At]3  [(Y^)At][(-^)At]2 


-[l-exp[-(YpL)At]]  [ 


[(Yp1’ )  At]  [( — — ) At  -  (Yp1 ) At] 2 

K  £o£r  k 

(— )At  , 

_ £o£r _ I 

[(Yr)At]2[(— )At-(Y^)At]2 

£o£r 


-expKYp™-)A>][l+-5i— ][ 


(— )At 

£o£r 


(Yp  >4.  lW»-)A,rt(i-)At-(Y*)4l]‘ 


(A.  113) 


At  T  T 


LX  l  it  5  11 

®  x  5  JdT  J dx’ J dx”  )(^ ) 2 exp[-(Yp L)(x -x" )]  exp[-(^-)(x-x’)l 


0  0  0 


=  (Atr 


-{I-—! 

.  1  4 


t(YNL)At][(_JL)At]  H  ( — — )At  [( — — )At)2  [( — — )At]3 


2[(-£i-)At] 

£o£r 


3(Y^)At  [(Y^)At][(-^)At]  [(Y^)At]4[(-^-)At-(Y^)At] 

M  £o£r  £o£r 

h[l-exp[-(-^-)At]]  [— - - + - - - 


[(— ^-)At]4  [(YoL)At]2[(— )At]2 


[(Y  p  L  )At] 2 1( — — )At  -  (Y^1 )  At] 


NL  \  A4.i2 


h[l-exp[-(Y^)At]][- 


2  [(— 2-)At] 

£o£r 


[(Y?L  )At]3[(-^)At-(Y^)At] 


NL  \ 


+  exp(-(Yp“-)A.][l+-si— ][ 


2  [(— 5-)At] 

Eo£r 


(Yp  )At  [(Y^)At]3[('^2-)At-(Y^L)AtJ 

K  EoEr  k 


(A.  114) 
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At 

(j£,o>x  s  JdT  exPHYp)t]exp[(-^-)x] 


=  At  exp[(-^-)At]  [ 


exp[-(Yp )  A  t]  -  exp[-( — — )At] 


£o£r 


£o£r 


(— 2-)At  -  (Yp  )At 


]■ 


(A.  115) 


£o£r 


(<i)x  s  fdx  f  dT’exp[-(Yp)(t-T’)]  exp[(-^-)x] 

J  J  K  £o£r 

0  0 

exp[(-^-)At]  f  1-expK— )At]  exp[-(Yp  )A  t]  -  exp[-(-^-)At]  - 


=  (At) 


Eo£r 


(Yp)At 


■{[■ 


_ £o£r 

(— )At 

£o£r 


£o£r 


(~*~)At  -  (Yp )  At 

£o£r 


]}.<A. 


116) 


At  T 


(Jtp2)x  s  fdxfdT’(^-)exp[-(Yp)(T-T,)]exp[(-^-)T] 
M  J  J  At  r  EoEr 


=  (At)2 


exp[(-^)At] 

£o£r 


{>-[ 


1 — exp[-( — — )At] 


£o£r 


[(■~^“)At][(Yp  )At]  (— )At 

£o£r  Eo£r 


][ 


i+ 


(— )A- 

£o£r 

(Yp)Af 


-[ 


exp[-(Yp )  A  t]  -  exp[-(-^)  At] r  (-^-)At . 


£o£r 


(— i_)At  -  (Yp  )At 


(Yp)At 


£o£r 


(A.  117) 


At 

(n£o)x  =  f  dx  exp[-(Y^L)x]exp[(-^-)x] 

«  y  £o£r 


:  At  exp[(-^-)At]  [■ 

EoEr 


exp[-(Y^ )  A  t]  -  exp[-(— ^-)At) 
_ _ EoEr 

(— )At-(Y^)At 
EoEr 


l 


(A.  11 8) 


At 

(Jt^-)x  =  I*  dx  x  exp[-(Y!fOx]exp[(— )(x)] 

J  K  £o£r 


=  (At)2 


exp[( — — )At] 

_ £o£r _ 

(£x_)At-(Y^)At 

£o£r  K 


exp[-(Y£IL  )A  t]  -exp[-( — — )At] 

exp[-(Y^  )A  t] - ^5. 

(— i-)At-(Y^fL)At 

£o£r  h 


]. 


(A.  119) 
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At  T 

(<k)x  =  fdx  fdT’exp[-(Y^L)(x-x,)]exp[(— )t] 

J  H  EoEr 


0  0 


exp[(— )At]  1-expK— )At]. 


=  (At) 


EoEr 


(YpNL)At 


EoEr 


(— )At 

EoEr 


-[ 


exp[-(Yp  L  )A  t]  -  exp[-(J1^-)At]  - 


EoEr 


(— 2-)At  -  (Y^  )At 

EoEr 


]}. 


At  T 


(<k>x  s  fdxfdx’(^)exp[-(Y^)(x-x’)]exp[(-^-)x] 
P  J  J  At  P  £o£r 


=  (At)2 


exp[(-^-)At] 

_ EoEr _ 

t(~^_)At][(Y^L  )At] 

EoEr  k 


{-[ 


1— exp[-< — — )At] 


EoEr 


(— )At 

EoEr 


][ 


i+ 


(-^)A- 

EoEr 


-[ 


exp[-(YpL)At]-exp[-(-^-)Atl  (— )At 

v  °  1  1  EoEr  I  l 


EoEr 


(-^)At-(Y^)At 

EoEr  ' 


'-][ 


(y?l)a f 


exp[(-^-)At]  f 

=  (At)2 - — - {  1 - - - 

[(~i_)At][(Y^rL  )At]  (— ^~)At 

EoEr  EoEr 


2 

(Y^)At 


exp[-( — — )At 

EoEr 


]][ 


[(-^-)At]2 

EoEr 


[(Y 


NL)At][(-^)At] 

EoEr 


f(YpL)At]2 


2  [(  — )At]  o  11 

-[ - — - J  Lexp[-(YpL  )A  t] - exp[-( — — ) At]  J  }, 

[(YpL)At]2t(— )Al-(Y^)At]  e°£R 


(A.  120) 


(A.  121) 


(A.  122) 
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At  T 


(7^5)  x  5  f  dx  f  dx’  (y- )  expHv^XT  -  x’)]  exp[(— f-)x] 

r’  J  J  At  £o£r 


exPt(  — )At]  f  r  r  1  i 

=  (At)2 - ^ - {  1  - [l ~ exp[-( — — )At] J  [— - J 

~  -  v  8o£r  n 


[(— 5-)At][(Y^IL  )At] 

8o8r 

(— )At 

.  I  8o£r 


(— )At 

8o8r 


'[(^)At-(Y^)At]2 

8o8r 


]  [expHY^  )A  t]  -  exp[-(— )At]] 


8o8r 


-[ 


8o£r  _ 


(^)At-(YpNL)At 

8o8r  f 


]  [exp[-(Y^1L)At]] 


(— )At 

8o8r 


EoEr 


[(ff-)At-(Y^)At]2 


-]  [exp[-(Y^  )A  t]  -  exp[-(— )At]] , 


8o8r 


(ti$)  x  =  Jdx  Jdx’ (■“)(“)  exp[-(Y^)(T  -  x’)]  exP[(^“)x] 

{ 


=  (At) 


exp[(-^)At] 

8o8r 


[(— 25-)At][(Y^IL  )At] 

8o8r 

-  [l  -  exp[-(-^-)  At]  ]  [ 

8o£r 


_i_ 

,NLs 


(-^)At  (Yp^)At 
EoEr 
2 


1 


[(— ^- )  At] 2  [(Yp  L )  At]  )At] 

8o8r  8o8r 


r  (  I  r  ~| 

+  [ - ^ - J  Lexp[-(Y^  )A  t]  J 

[(Y^)At][(— )At-(Y^)At] 

P  8o8r 

-[ - — - ]  [exp[-(Y^  )A  t]  -exp[-(-^-)At]]  } , 

[(Y^)At][(-^)At-(Y^)At]2  £°eR 


(A.  123) 


(A.  124) 
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2 


At  T  ? 

Ox  *  I’d*  f dx’(-f  )(^7 )2exp[-(YpNL)(x-T’)]  exp[(-^-)x] 
J  J  At  At  K  £o£r 


:(At)2 


At  At 

o 

exp[(-^*-)At] 

_ £o£r _ 

0)At][0)At] 

£o£r  k 


{- 


(— )At  [(— *-)At] 

£o£r  £o£r 


2  (Y^)At 


[( — ~)At]2  t(Yp  )At][( — — )At] 


£o£r 


Eo£r 


-[l-exp[-(-^-)At]]  [— - - +- 


£o£r 


[( — — )At]3  [(Y  p  ) At]  [( — — )At] 1  [(Y  p  )At]  [( — — )At] 

£o£r  £o£r  £o£r 


2  [(— 2t—)At] 

£o£r 


"]  [exp[_(YpL )A  t]] 


[(Yp"L)At][(^-)At-(YpNL)At] 

£o£r 

2  )At] 

£o£r 


[(YpL)At][(— i-)At  -  (y^  )At]2 

K  £o£r  p 


■]  [exp[-(Y^L)At]-exp[-(-^-)At]]  },  (A.125) 


£o£r 


For  the  v  component: 

Replace  x—»y  of  matrix  elements  defined  for  the  x  component  above. 
For  the  z  component: 

Replace  x^z  of  matrix  elements  defined  for  the  x  component  above. 
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Flow  Chart 
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